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We report on numerical results from an independent for- 
malism to describe the quasi-equilibrium structure of nonsyn- 
chronous binary neutron stars in general relativity. This is an 
important independent test of controversial numerical hydro- 
dynamic simulations which suggested that nonsynchronous 
neutron stars in a close binary can experience compression 
prior to the last stable circular orbit. We show that, for com- 
pact enough stars the interior density increases slightly as 
irrotational binary neutron stars approach their last orbits. 
The magnitude of the effect, however, is much smaller than 
that reported in previous hydrodynamic simulations. 

PACS number(s): 97.80.Fk, 04.25.Dm, 04.40.Dg, 97.60. Jd 



I. INTRODUCTION 

The physical processes occurring during the last orbits 
of a neutron-star binary are a subject of intense current 
debate |l|-[rj]]. In part, this recent surge in interest stems 
from relativistic numerical hydrodynamic simulations in 
which it has been noted |0-S| that as the stars approach 
each other their interior density increases. Indeed, for an 
appropriate mass and equation of state, previous numeri- 
cal simulations indicated that binary neutron stars might 
collapse individually toward black holes seconds prior to 
merger. This compression effect would have a significant 
impact on the anticipated gravity- wave signal from merg- 
ing neutron stars and may also provide an energy source 
for cosmological gamma-ray bursts . 

In view of the unexpected nature of this neutron star 
compression effect, and its possible repercussions, as well 
as the extreme complexity of strong field general rela- 
tivistic hydrodynamics, it is of course imperative that 
there be an independent confirmation of the existence of 
neutron star compression before one can be convinced of 
its operation in binary systems. This is particularly im- 
portant since many authors p JlOfl have searched for but 
not observed this effect in Newtonian tidal forces [[| , first 
post-Newtonian (1PN) dynamics |3 10 1, tidal expansions 
|0,gj, or in binaries in which rigid corotation has been 
imposed |^,[l^] . In ref . Q] it has been argued that none of 
the above works could or should have observed the effect, 
since the compression only dominates over tidal forces for 
stars with a realistic compaction ratio [(Mg/R)oo\ that 
are in a nearly irrotational hydrodynamic state (i.e. little 



spin relative to a distant observer). Indeed, irrotational 
stars may be a likely configuration near the final orbits 
as corotation would demand an unrealistically large vis- 
cosity in neutron stars |L3j . In hydrodynamic relaxation 
calculations || the stars even seem to prefer a nearly ir- 
rotational state. Rasio and Shapiro JL4| provide a recent 
review of the state of the aforementioned controversy. 

With this in mind it is of particular interest that a new 
formalism [|15| |19[ has been proposed in which the hydro- 
static quasi-equilibrium of irrotational stars in a binary 
can be solved independently of the complexities of (3+1) 
nunrerical relativistic hydrodynamics. Thus, this pro- 
vides an opportunity to independently test the hydrody- 
namics result. Previous results using this formalism have 
been presented by Bonazzola et al. |1T[ . They show that 
for an irrotational binary system composed of neutron 
stars with {Mg/R)oo of 0.14 there is almost no evidence 
of the compression effect (although they did note that the 
central density for irrotational stars remains much higher 
than for corotating stars as the orbit decays). In this pa- 
per we report results for two different systems involving 
stars with {Mq/R)^ of 0.14 and 0.19. We show that, 
while the former sequence shows almost no compression 
in agreement with pdj |, stars with a higher compaction 
ratio seem to experience a slight central compression as 
the stars approach. 

We also note that these sequences have an important 
intrinsic value beyond the controversy about the com- 
pression effect. They provide realistic solutions to the ini- 
tial value problem for neutron-star binary systems that 
can be used as starting points of fully dynamical rela- 
tivistic hydrodynamical simulations. The fact that they 
are valid in the strong field regime makes them more 
attractive than post-Newtonian counterparts, since the 
simulations can be initiated at stages very close to the 
final merger of the stars. 



II. THE METHOD 

The method we use to determine the internal structure 
of stars in irrotational quasi-equilibrium configurations is 
essentially that originally proposed by Bonazzola, Gour- 
goulhon & Marck lla] and as simplified by Teukolsky pw . 
The derivation of the relevant equations can be found in 
those papers. The essential ingredient of this approach is 



that, if the fluid vorticity is zero (e.g. as in irrotational 
stars) , the specific momentum density per baryon can be 
written as the gradient of a scalar potential, 



hUn = Vnt/J 



(1) 



where u^ is the covariant four velocity and h is the rela- 
tivistic enthalpy, h = 1 + e + — , where e is the internal 
energy per unit of baryon rest mass, P is the pressure, 
and po = rnBTiB is the proper baryon rest mass density, 
with tib the baryon number density. The potential ij> can 
be obtained from the solution of a Poisson-likc equation: 



D*D^j = D^ 



°'*->h<i^r 



(2) 



where Di are spatial covariant derivatives, and 

X = C + B : >D^^a[h 2 + {D l ^) 2 } 1 / 2 , (3) 

where C is a constant. The quantities a and B l are 
obtained from the usual ADM (3+1) metric 



ds 1 



-(a 2 - f3i(3 l )dt 2 + 2(3idx l dt + j^dx'dx 3 , (4) 



such the B % is the shift vector in the rotating frame, B l = 
ff + (u> x r) 1 , where u> is the angular velocity of orbital 
motion. 

Equation (0) must be solved by imposing a boundary 
condition at the stellar surface, 



D % i/>- 



X 



B l }D z n B 



surf 



= o 



(5) 



For irrotational stars the Bernoulli integral for the matter 
distribution then becomes: 



fc 2 = -(I>»ZW+4 
Or 



(6) 



Equation (ra) uniquely determines the equilibrium struc- 
ture of the stars. 

We also compute stars in constant corotation. In this 
case the relativistic Bernoulli equation can be written 
|,H h/u° = constant. 

In the present work we consider a polytropic equation 
of state, P = Kpl, with r = 2 and K = 1.13 x 10 5 erg 
cm 3 g -2 . This gives a maximum neutron-star gravita- 
tional mass of 1.46 M . In these simulations we con- 
sider two equal-mass neutron stars in two different se- 
quences of stable orbits. The first corresponds to stars 
with baryon mass of Mb = 1-29 Mq, gravitational mass 
in isolation of Mqoc = 1-21 Mq, and compaction ratio 
of (Mg/R)oc = 0.14. This sequence is comparable to 
the one studied by Bonazzola et al. ill]. The second se- 
quence corresponds to a systems of more compact stars, 
with baryon mass of Mb — 1-55 Mq, gravitational mass 



in isolation of Mqoo — 1.41 Mq, and compaction ratio 
of 0.19. For the grid resolution of this study (~ 40 zones 
across the star) we obtain a proper central density in 
isolation of p x = 0.96 x 10 14 g cm -3 for the first se- 
quence and poo = 1.83 x 10 14 g cm~ 3 for the second one. 
As pointed out in it is important to study realistically 
compact neutron stars. Otherwise Newtonian tidal forces 
can dominate over the relativistic effects one desires to 
probe. 

The Einstein field equations are solved by imposing 
a conformally flat condition (CFC) on the three metric 
plpl. That is, the spatial three metric is constrained to 
be represented by a position dependent conformal factor 
<^> 4 times the Kronecker delta, 7y = (jfiSij . This is a com- 
mon metric choice for solving the initial value problem in 
numerical relativity. It is consistent with the approxima- 
tion of quasi-equilibrium circular orbits |4|Jl6|| which we 
wish to evaluate (see ref. Q for a detailed discussion). 

The advantage of this method is that the determina- 
tion of the metric coefficients reduces to the solution of 
flat-space Poisson like equations M. For example, using 
the Hamiltonian constraint EG] in combination with the 
maximal slicing condition tr(K) = 0, the equation for <p 
becomes 



V z 6 = -An— 
2 



Po hW 2 - P + — KijK 1 ' 

167T 



(7) 



where, W = au°. A similar equation can be written |2j 
for the lapse function. 

The shift vector (3 l is further decomposed [pi[ : f3 l = 
G l — (1/4) V J x, and introduced into the ADM momentum 
constraint equation to obtain two elliptic equations 



V 2 X = V t G l 



V 2 G i = 2V ln(a(j>- & )K ij - Uira^S 1 



(8) 



(9) 



These equations have been solved with boundary condi- 
tions provided by the first terms in a multipole expan- 
sion of the fields as described in |]12). An equation for 
the extrinsic curvature K l i follows |2| from the ADM 
momentum constraint equation with the maximal slicing 
condition. 

This set of equations (g), (||)-(||) is solved numerically 
using an iterative algorithm based upon a specially de- 
signed elliptic solver. This method consists of a combina- 
tion of multigrid algorithms and domain decomposition 
techniques |l2|]23J l and utilizes a code which was devel- 
oped independently of the hydrodynamics code of |l]||] . 

Solutions are obtained for specific values of the coordi- 
nate distance between stars and the total baryonic mass. 
In this way we can construct a constant baryonic-mass 
sequence of orbits with a minimum number of code runs. 
This sequence is a collection of semi-stable orbits that 
are connected by the inspiral motion of the stars. 



The problem consists essentially in the numerical solu- 
tion of a set of elliptic equations. The cases of corotating 
and irrotational systems share the same set of equations 
for the metric fields. The irrotational systems, however, 
demand the solution of an extra elliptic equation (0) for 
the description of the stellar internal structure. This 
poses a very special problem since the boundary con- 
dition (||) is to be satisfied on the stellar surface and not 
on the grid boundaries as for the rest of the elliptic equa- 
tions. This is challenging to implement numerically on a 
fixed Eulerian grid, since the stellar surface is a spheroid 
embedded in a Cartesian grid. We refer the reader to J2JJ 
for the details of our approach to this calculation. Note, 
however, that this is a completely different numerical ap- 
proach to that of Bonazzola et al. who utilized spectral 
methods to solve the elliptic equations. 
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FIG. 1. Change in central density relative to a single iso- 
lated star (Ap/p) as a function of the coordinate distance. 
The black (white) circles represent our calculations for se- 
quences with (Ma/R)oo = 0.19 (0.14). The squares show the 
result obtained by Bonazzola et al. for a sequence composed 
of stars with (M G /R)oo = 0.14. 

To quantify the numerical accuracy, the code was 
tested against the Newtonian irrotational sequences ob- 
tained by Uryu and Eriguchi p4j in two different ways. 
In one test, the code was stripped of all relativistic terms 
to reduce it to Newtonian physics. In a second test, or- 
bits for low mass stars (Mb — 0.10 M Q ) were calcu- 
lated using the relativistic code to approach the Newto- 
nian regime (Mq/R ~ 0.01). Comparisons were done 
for orbits with fixed separation distances between stel- 
lar centers (d in the units of |24] ). Most importantly, 
both of these tests exhibit the Newtonian expectation 
of no discernible change in central density relative to 
that of isolated stars, with (Ap/p <0.1%). The result- 



and fl = 0.241 (3%); In the second test, d = 3.88, and 
E = -1.082 (4%), J = 1.413 (1%), and f2 = 0.203 (6%). 
These results are consistent with the expected numerical 
accuracy of this comparison. The code was also tested 
simulating corotating binary orbits (see |l2[). 

Finally, Flanagan 122] pointed out an inconsistency in 
the definition of momentum density in the original hydro- 
dynamics calculations [0-pl ■ This problem is not present 
in these calculations. 



III. RESULTS 

For the present study we have found solutions to the 
initial value equations described above for semi-stable cir- 
cular orbits for a binary system of identical neutron stars 
extending from the post-Newtonian regime to the inner- 
most orbit for which we obtain a stable solution. 

Figure [j] shows the fractional change Ap/p in the 
proper central rest-mass density relative to the central 
density p^ of an isolated star of the same resolution. 
Results are plotted as a function of the normalized coor- 
dinate distance; i. e. d/M 1A = d/(M G /\A M ). 
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ing values (and percentage difference from results of |24|) 
for total energy, total angular momentum, and orbital 
angular frequency were as follows: For the first test, 
d = 3.53, and E = -1.142 (0.1%), J = 1.278 (4%), 
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FIG. 2. Binding energy defined as (M G - M Goo )/M Goo , 
where Mq is half the ADM mass of the system, vs coordinate 
distance. The black (white) circles correspond to the sequence 
of stable orbits with (M G /R)oo = 0.19 (0.14). The squares 
show the result obtained by Bonazzola et al. for a sequence 
composed of stars with (M G /R)oo = 0.14. 

The squares show the orbits obtained by Bonazzola et 
al. |ll]] for a system composed of stars with (Mg/R)<x> — 
0.14. The black (white) circles represent our calculations 
for sequences with (Mq/R)oc = 0.19 (0.14). The numer- 
ical error in Ap/p for our irrotational points is conser- 
vatively estimated to be ~ ±0.005 based upon the code 
tests here described as well as various convergence tests. 

Note, that both the Bonazzola et al. results and ours 
for (Mg/R)oo — 0.14 show only a very slight increase 



in the central density before the tidal effects begin to 
dominate at short separations. The difference between 
the two curves is probably attributable to differences in 
the numerical algorithms applied (i.e. coordinate mesh, 
boundary conditions, elliptic solver, etc). Nevertheless 
they are consistent to within numerical accuracy. We em- 
phasize that our point is not about quantitative details 
at the level of ±0.005 in Ap/p, but about the qualitative 
trend of increasing central density as the compaction ra- 
tio increases. 

The sequence with (Mg/R)oo = 0.19, clearly, raises 
above the zero line, ending in a maximum value of ~ 1.5% 
increase in the central density for the last stable orbit 
found. This result agrees with a similar trend of increas- 
ing compression with increasing compaction ratio, pre- 
sented recently by Bonazzola et al. |25). This trend de- 
rives from the fact that the compression effect competes 
with the tidal deformation and the latter is stronger for 
more extended (i.e. smaller compaction ratio) stars. 

In both cases the central density approaches the 
isolated-star limit at large orbital separations. We have 
checked that the same trend emerges in a plot of the av- 
erage density, so this effect could not be an artifact of 
the stellar center being a special point. 

As a final point, figure @ shows the relative binding en- 
ergy of the system defined as {Mq — Mgoc)/Mgoo, where 
Mq is half the ADM mass of the system. We note that 
the relative binding energy does not depend strongly on 
the compaction ratio within the numerical accuracy of 
our results. Another feature of these curves is that they 
lack a turning point in the sequence. This is in qualitative 
agreement with the Newtonian irrotational sequences for 
polytropic index n — 1 obtained by Uryu and Eriguchi 

In summary, the present independent study has ob- 
tained the qualitative result of increasing central density 
as an irrotational binary orbit decays, for neutron stars 
with compaction ratio of 0.19. The magnitude of the ef- 
fect, however, is significantly less than that of the previ- 
ously reported hydrodynamic results [|l|-|3| . Preliminary 
recent results [|26| indicate that such reduction in com- 
pression effect is consistent with the effect of a correction 
proposed by Flanagan |2^] when applied to the hydrody- 
namic simulations. 

Although it is evident that some compression effect ex- 
ists, it is not yet clear whether this remaining effect is real 
or a consequence of the conformally flat metric approxi- 
mation. The magnitude of the apparent effect noted here 
is roughly comparable to the possible uncertainty intro- 
duced by this approximation (cf. 0). The ultimate test 
of whether the effect is real will therefore require an accu- 
rate fully dynamical relativistic treatment. Such a test is 
hopefully coming from the neutron-star Grand Challenge 
collaboration. 
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